﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Fortuitus.Core
{
    public static class Statistics
    {
        /// <summary>
        /// Calculates the fold change of two number greater than zero.
        /// </summary>
        /// <param name="a">The first number.</param>
        /// <param name="b">The second number.</param>
        /// <returns>The fold change that gives the number of time b is greater than a.</returns>
        /// <exception cref="ArgumentException">Both numbers need to be greater than zero.</exception>
        public static double FoldChange(double a, double b)
        {
            if (a <= 0 || b <= 0)
                throw new ArgumentException("Foldchange is only defined for number greater than zero");

            double foldchange;
            double ratio = b / a;
            if (ratio < 1)
                foldchange = 1 / ratio;
            else
                foldchange = ratio;

            return foldchange;

        }

        public static double Welchs_t_test(IVector a, IVector b, out int  df)
        {
            double n1 = a.Length;
            double n2 = b.Length;
            double x1 = a.Mean();
            double x2 = b.Mean();
            double s1 = a.VarianceS();
            double s2 = b.VarianceS();
           
            double Sx1x2 = Math.Sqrt((s1/n1)+(s2/n2));
            double t = (x1-x2)/Sx1x2;
            df = WelchSatterthwaiteEquation(s1, s2, n1, n2);

            return t;

        }

        public static double t_test(IVector a, IVector b, out int df)
        {
            double n1 = a.Length;
            double n2 = b.Length;
            double x1 = a.Mean();
            double x2 = b.Mean();
            double s1 = a.VarianceS();
            double s2 = b.VarianceS();
            df = (int)n1 + (int)n2 - 2;
            double Sx1x2 = Math.Sqrt( (((n1 - 1) * s1) + ((n2 - 1) * s2))/df );
            double t = (x1 - x2) / (Sx1x2 * Math.Sqrt((1/n1)+(1/n2))) ;
           

            return t;

        }


        public static double OneWayANOVA(IEnumerable<IVector> vectors, out int df1, out int df2)
        {
            var vectorList = vectors.ToList();
            var overallMean = (double)vectorList.Sum(x => x.Sum(y => y)) / (double)vectorList.Sum(x => x.Length);
            double SSB = 0; //SSB sum of squares between groups
            for (int i = 0; i < vectorList.Count; i++)
            {
                double mean = vectorList[i].Mean();
                double n = vectorList[i].Length;
                SSB += (n * Math.Pow((mean - overallMean),2));
            }

            df1 = vectorList.Count - 1;
            var MSB = SSB / df1; //MSB mean square between groups

            double SSW = 0; // SSW sum of squares within groups
            for (int i = 0; i < vectorList.Count; i++)
            {
                var vector = vectorList[i].ToEditableVector();
                double mean = vector.Mean();
                //center
                vector.Apply(x => x - mean);
                //square
                vector.Apply(x => Math.Pow(x,2));
                SSW += vector.Sum();
            }
            df2 = vectorList.Sum(x => x.Length) - vectorList.Count;
            double MSW = SSW / df2; //MSW mean square within groups
            double F = MSB / MSW;
            return F;

        }

        public static int WelchSatterthwaiteEquation(double s1, double s2, double n1, double n2)
        {
            double A = ((s1 / n1) + (s2 / n2)) * ((s1 / n1) + (s2 / n2));
            double B = (((s1 / n1) * (s1 / n1)) / (n1 - 1)) + (((s2 / n2) * (s2 / n2)) / (n2 - 1));
            double ret = A / B;
            return (int)Math.Floor(ret);
        }

        public static double Students_PvalueTwoTailed(double tvalue, int df)
        {
            if (double.IsNaN(tvalue))
                return double.NaN;
            return (IBeta(((double)df / 2), 0.5, ((double)df / ((double)df + Math.Pow(tvalue, 2)))));
        }

        public static double Students_PvalueOneTailed
            (double tvalue, int df)
        {
            return Students_PvalueTwoTailed(tvalue,df)/2;
        }

        public static double F_Pvalue(double fvalue, int df1, int df2)
        {
            if (double.IsNaN(fvalue))
                return double.NaN;

            double x = ((df1 * fvalue)/df2) / (1+ ( (df1 * fvalue)/df2));

            return 1 - IBeta((double)df1 / 2, (double)df2 / 2, x);
        }

        public static int Factorial(int number)
        {
            if (number == 0 || number == 1 )
                return 1;
            else
                return number * (Factorial(number - 1));

        }

        /// <summary>
        /// Calculates the Pearson Correlation Coefficient between to vectors.
        /// </summary>
        /// <param name="a">The first vector.</param>
        /// <param name="b">The second vector.</param>
        /// <returns>The Correlation Coefficient.</returns>
        public static double PearsonCorrelationCoefficient(IVector a, IVector b)
        {
            if (a.Count() != b.Count())
                throw new ArgumentException("Vector a has to be of the same size as vector b");
            if(a.IsEmptyOrHasNoValues() || b.IsEmptyOrHasNoValues())
                throw new ArgumentException("Vectors can't be empty or contain just missing values");

            int n = a.Count();
            double meanA = a.Mean();
            double meanB = b.Mean();
            double stdDevA = a.StandardDevP();
            double stdDevB = b.StandardDevP();

            double pc = 0;
            for (int i = 0; i < n; i++)
            {
                double A = (a[i] - meanA) / stdDevA;
                double B = (b[i] - meanB) / stdDevB;
                pc += A * B;
            }
            
            pc = pc / (n);
            return pc;
        }

        /// <summary>
        /// Calculates the Spearman Correlation Coefficient between to vectors.
        /// </summary>
        /// <param name="a">The first vector.</param>
        /// <param name="b">The second vector.</param>
        /// <returns>The Correlation Coefficient.</returns>
        public static double SpearmanCorrelationCoefficient(IVector a, IVector b)
        {
            if (a.Count() != b.Count())
                throw new ArgumentException("Vector a has to be of the same size as vector b");
            if (a.IsEmptyOrHasNoValues() || b.IsEmptyOrHasNoValues())
                throw new ArgumentException("Vectors can't be empty or contain just missing values");

            int n = a.Count();
            IVector ranksA = a.ToRanks();
            IVector ranksB = b.ToRanks();

            //var test = PearsonCorrelationCoefficient(ranksA, ranksB);
            //return test;
            double meanA = ranksA.Mean();
            double meanB = ranksB.Mean();

            //rho = 
            //           X
            // ----------------------------
            //     sqrt(Y1 * Y2)

            double X = 0;
            for (int i = 0; i < n; i++)
            {
                double A = (ranksA[i] - meanA);
                double B = (ranksB[i] - meanB);

                X += A * B;
            }
            double Y1 = 0;
            for (int i = 0; i < n; i++)
            {
                double A = (ranksA[i] - meanA);
                Y1 +=  A * A;
            }

            double Y2 = 0;
            for (int i = 0; i < n; i++)
            {
                double A = (ranksB[i] - meanB);
                Y2 += A * A;
            }

            var Y = Math.Sqrt(Y1 * Y2);

            var spc = X / Y;

            return spc;
        }
    

        /// <summary>
        /// Returns the incomplete beta function evaluated from zero to xx.
        /// </summary>
        /// <param name="aa"></param>
        /// <param name="bb"></param>
        /// <param name="xx"></param>
        /// <returns></returns>
        public static double IBeta(double aa, double bb, double xx)
        {
            double a, b, t, x, xc, w, y;
            bool flag;

            if (aa <= 0.0 || bb <= 0.0) throw new
                                            ArithmeticException("ibeta: Domain error!");

            if ((xx <= 0.0) || (xx >= 1.0))
            {
                if (xx == 0.0) return 0.0;
                if (xx == 1.0) return 1.0;
                throw new ArithmeticException("ibeta: Domain error!");
            }

            flag = false;
            if ((bb * xx) <= 1.0 && xx <= 0.95)
            {
                t = pseries(aa, bb, xx);
                return t;
            }

            w = 1.0 - xx;

            /* Reverse a and b if x is greater than the mean. */
            if (xx > (aa / (aa + bb)))
            {
                flag = true;
                a = bb;
                b = aa;
                xc = xx;
                x = w;
            }
            else
            {
                a = aa;
                b = bb;
                xc = w;
                x = xx;
            }

            if (flag && (b * x) <= 1.0 && x <= 0.95)
            {
                t = pseries(a, b, x);
                if (t <= MACHEP) t = 1.0 - MACHEP;
                else t = 1.0 - t;
                return t;
            }

            /* Choose expansion for better convergence. */
            y = x * (a + b - 2.0) - (a - 1.0);
            if (y < 0.0)
                w = incbcf(a, b, x);
            else
                w = incbd(a, b, x) / xc;

            /* Multiply w by the factor
                   a      b   _             _     _
                  x  (1-x)   | (a+b) / ( a | (a) | (b) ) .   */

            y = a * Math.Log(x);
            t = b * Math.Log(xc);
            if ((a + b) < MAXGAM && Math.Abs(y) < MAXLOG && Math.Abs(t) < MAXLOG)
            {
                t = Math.Pow(xc, b);
                t *= Math.Pow(x, a);
                t /= a;
                t *= w;
                t *= Gamma(a + b) / (Gamma(a) * Gamma(b));
                if (flag)
                {
                    if (t <= MACHEP) t = 1.0 - MACHEP;
                    else t = 1.0 - t;
                }
                return t;
            }
            /* Resort to logarithms.  */
            y += t + LGamma(a + b) - LGamma(a) - LGamma(b);
            y += Math.Log(w / a);
            if (y < MINLOG)
                t = 0.0;
            else
                t = Math.Exp(y);

            if (flag)
            {
                if (t <= MACHEP) t = 1.0 - MACHEP;
                else t = 1.0 - t;
            }
            return t;
        }

        /// <summary>
        /// Returns the continued fraction expansion #1 for incomplete beta integral.
        /// </summary>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="x"></param>
        /// <returns></returns>
        private static double incbcf(double a, double b, double x)
        {
            double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
            double k1, k2, k3, k4, k5, k6, k7, k8;
            double r, t, ans, thresh;
            int n;
            double big = 4.503599627370496e15;
            double biginv = 2.22044604925031308085e-16;

            k1 = a;
            k2 = a + b;
            k3 = a;
            k4 = a + 1.0;
            k5 = 1.0;
            k6 = b - 1.0;
            k7 = k4;
            k8 = a + 2.0;

            pkm2 = 0.0;
            qkm2 = 1.0;
            pkm1 = 1.0;
            qkm1 = 1.0;
            ans = 1.0;
            r = 1.0;
            n = 0;
            thresh = 3.0 * MACHEP;
            do
            {
                xk = -(x * k1 * k2) / (k3 * k4);
                pk = pkm1 + pkm2 * xk;
                qk = qkm1 + qkm2 * xk;
                pkm2 = pkm1;
                pkm1 = pk;
                qkm2 = qkm1;
                qkm1 = qk;

                xk = (x * k5 * k6) / (k7 * k8);
                pk = pkm1 + pkm2 * xk;
                qk = qkm1 + qkm2 * xk;
                pkm2 = pkm1;
                pkm1 = pk;
                qkm2 = qkm1;
                qkm1 = qk;

                if (qk != 0) r = pk / qk;
                if (r != 0)
                {
                    t = Math.Abs((ans - r) / r);
                    ans = r;
                }
                else
                    t = 1.0;

                if (t < thresh) return ans;

                k1 += 1.0;
                k2 += 1.0;
                k3 += 2.0;
                k4 += 2.0;
                k5 += 1.0;
                k6 -= 1.0;
                k7 += 2.0;
                k8 += 2.0;

                if ((Math.Abs(qk) + Math.Abs(pk)) > big)
                {
                    pkm2 *= biginv;
                    pkm1 *= biginv;
                    qkm2 *= biginv;
                    qkm1 *= biginv;
                }
                if ((Math.Abs(qk) < biginv) || (Math.Abs(pk) < biginv))
                {
                    pkm2 *= big;
                    pkm1 *= big;
                    qkm2 *= big;
                    qkm1 *= big;
                }
            } while (++n < 300);

            return ans;
        }


        /// <summary>
        /// Returns the continued fraction expansion #2 for incomplete beta integral.
        /// </summary>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="x"></param>
        /// <returns></returns>
        private static double incbd(double a, double b, double x)
        {
            double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
            double k1, k2, k3, k4, k5, k6, k7, k8;
            double r, t, ans, z, thresh;
            int n;
            double big = 4.503599627370496e15;
            double biginv = 2.22044604925031308085e-16;

            k1 = a;
            k2 = b - 1.0;
            k3 = a;
            k4 = a + 1.0;
            k5 = 1.0;
            k6 = a + b;
            k7 = a + 1.0;
            ;
            k8 = a + 2.0;

            pkm2 = 0.0;
            qkm2 = 1.0;
            pkm1 = 1.0;
            qkm1 = 1.0;
            z = x / (1.0 - x);
            ans = 1.0;
            r = 1.0;
            n = 0;
            thresh = 3.0 * MACHEP;
            do
            {
                xk = -(z * k1 * k2) / (k3 * k4);
                pk = pkm1 + pkm2 * xk;
                qk = qkm1 + qkm2 * xk;
                pkm2 = pkm1;
                pkm1 = pk;
                qkm2 = qkm1;
                qkm1 = qk;

                xk = (z * k5 * k6) / (k7 * k8);
                pk = pkm1 + pkm2 * xk;
                qk = qkm1 + qkm2 * xk;
                pkm2 = pkm1;
                pkm1 = pk;
                qkm2 = qkm1;
                qkm1 = qk;

                if (qk != 0) r = pk / qk;
                if (r != 0)
                {
                    t = Math.Abs((ans - r) / r);
                    ans = r;
                }
                else
                    t = 1.0;

                if (t < thresh) return ans;

                k1 += 1.0;
                k2 -= 1.0;
                k3 += 2.0;
                k4 += 2.0;
                k5 += 1.0;
                k6 += 1.0;
                k7 += 2.0;
                k8 += 2.0;

                if ((Math.Abs(qk) + Math.Abs(pk)) > big)
                {
                    pkm2 *= biginv;
                    pkm1 *= biginv;
                    qkm2 *= biginv;
                    qkm1 *= biginv;
                }
                if ((Math.Abs(qk) < biginv) || (Math.Abs(pk) < biginv))
                {
                    pkm2 *= big;
                    pkm1 *= big;
                    qkm2 *= big;
                    qkm1 *= big;
                }
            } while (++n < 300);

            return ans;
        }

        /// <summary>
        /// Returns the natural logarithm of gamma function.
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        public static double LGamma(double x)
        {
            double p, q, w, z;

            double[] A = {
						 8.11614167470508450300E-4,
						 -5.95061904284301438324E-4,
						 7.93650340457716943945E-4,
						 -2.77777777730099687205E-3,
						 8.33333333333331927722E-2
					 };
            double[] B = {
						 -1.37825152569120859100E3,
						 -3.88016315134637840924E4,
						 -3.31612992738871184744E5,
						 -1.16237097492762307383E6,
						 -1.72173700820839662146E6,
						 -8.53555664245765465627E5
					 };
            double[] C = {
						 /* 1.00000000000000000000E0, */
						 -3.51815701436523470549E2,
						 -1.70642106651881159223E4,
						 -2.20528590553854454839E5,
						 -1.13933444367982507207E6,
						 -2.53252307177582951285E6,
						 -2.01889141433532773231E6
					 };

            if (x < -34.0)
            {
                q = -x;
                w = LGamma(q);
                p = Math.Floor(q);
                if (p == q) throw new ArithmeticException("lgam: Overflow");
                z = q - p;
                if (z > 0.5)
                {
                    p += 1.0;
                    z = p - q;
                }
                z = q * Math.Sin(Math.PI * z);
                if (z == 0.0) throw new
                                  ArithmeticException("lgamma: Overflow");
                z = LOGPI - Math.Log(z) - w;
                return z;
            }

            if (x < 13.0)
            {
                z = 1.0;
                while (x >= 3.0)
                {
                    x -= 1.0;
                    z *= x;
                }
                while (x < 2.0)
                {
                    if (x == 0.0) throw new
                                      ArithmeticException("lgamma: Overflow");
                    z /= x;
                    x += 1.0;
                }
                if (z < 0.0) z = -z;
                if (x == 2.0) return Math.Log(z);
                x -= 2.0;
                p = x * polevl(x, B, 5) / p1evl(x, C, 6);
                return (Math.Log(z) + p);
            }

            if (x > 2.556348e305) throw new
                                      ArithmeticException("lgamma: Overflow");

            q = (x - 0.5) * Math.Log(x) - x + 0.91893853320467274178;
            if (x > 1.0e8) return (q);

            p = 1.0 / (x * x);
            if (x >= 1000.0)
                q += ((7.9365079365079365079365e-4 * p
                    - 2.7777777777777777777778e-3) * p
                    + 0.0833333333333333333333) / x;
            else
                q += polevl(p, A, 4) / x;
            return q;
        }

        /// <summary>
        /// Evaluates polynomial of degree N with assumtion that coef[N] = 1.0
        /// </summary>
        /// <param name="x"></param>
        /// <param name="coef"></param>
        /// <param name="N"></param>
        /// <returns></returns>		
        private static double p1evl(double x, double[] coef, int N)
        {
            double ans;

            ans = x + coef[0];

            for (int i = 1; i < N; i++)
            {
                ans = ans * x + coef[i];
            }

            return ans;
        }

        /// <summary>
        /// Returns the power series for incomplete beta integral. Use when b*x is small and x not too close to 1.
        /// </summary>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="x"></param>
        /// <returns></returns>
        private static double pseries(double a, double b, double x)
        {
            double s, t, u, v, n, t1, z, ai;

            ai = 1.0 / a;
            u = (1.0 - b) * x;
            v = u / (a + 1.0);
            t1 = v;
            t = u;
            n = 2.0;
            s = 0.0;
            z = MACHEP * ai;
            while (Math.Abs(v) > z)
            {
                u = (n - b) * x / n;
                t *= u;
                v = t / (a + n);
                s += v;
                n += 1.0;
            }
            s += t1;
            s += ai;

            u = a * Math.Log(x);
            if ((a + b) < MAXGAM && Math.Abs(u) < MAXLOG)
            {
                t = Gamma(a + b) / (Gamma(a) * Gamma(b));
                s = s * t * Math.Pow(x, a);
            }
            else
            {
                t = LGamma(a + b) - LGamma(a) - LGamma(b) + u + Math.Log(s);
                if (t < MINLOG) s = 0.0;
                else s = Math.Exp(t);
            }
            return s;
        }


        /// <summary>
        /// Returns the gamma function of the specified number.
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        public static double Gamma(double x)
        {
            double[] P = {
						 1.60119522476751861407E-4,
						 1.19135147006586384913E-3,
						 1.04213797561761569935E-2,
						 4.76367800457137231464E-2,
						 2.07448227648435975150E-1,
						 4.94214826801497100753E-1,
						 9.99999999999999996796E-1
					 };
            double[] Q = {
						 -2.31581873324120129819E-5,
						 5.39605580493303397842E-4,
						 -4.45641913851797240494E-3,
						 1.18139785222060435552E-2,
						 3.58236398605498653373E-2,
						 -2.34591795718243348568E-1,
						 7.14304917030273074085E-2,
						 1.00000000000000000320E0
					 };

            double p, z;

            double q = Math.Abs(x);

            if (q > 33.0)
            {
                if (x < 0.0)
                {
                    p = Math.Floor(q);
                    if (p == q) throw new ArithmeticException("gamma: overflow");
                    //int i = (int)p;
                    z = q - p;
                    if (z > 0.5)
                    {
                        p += 1.0;
                        z = q - p;
                    }
                    z = q * Math.Sin(Math.PI * z);
                    if (z == 0.0) throw new ArithmeticException("gamma: overflow");
                    z = Math.Abs(z);
                    z = Math.PI / (z * stirf(q));

                    return -z;
                }
                else
                {
                    return stirf(x);
                }
            }

            z = 1.0;
            while (x >= 3.0)
            {
                x -= 1.0;
                z *= x;
            }

            while (x < 0.0)
            {
                if (x == 0.0)
                {
                    throw new ArithmeticException("gamma: singular");
                }
                else if (x > -1.0E-9)
                {
                    return (z / ((1.0 + 0.5772156649015329 * x) * x));
                }
                z /= x;
                x += 1.0;
            }

            while (x < 2.0)
            {
                if (x == 0.0)
                {
                    throw new ArithmeticException("gamma: singular");
                }
                else if (x < 1.0E-9)
                {
                    return (z / ((1.0 + 0.5772156649015329 * x) * x));
                }
                z /= x;
                x += 1.0;
            }

            if ((x == 2.0) || (x == 3.0)) return z;

            x -= 2.0;
            p = polevl(x, P, 6);
            q = polevl(x, Q, 7);
            return z * p / q;

        }


        /// <summary>
        /// Return the gamma function computed by Stirling's formula.
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        private static double stirf(double x)
        {
            double[] STIR = {
							7.87311395793093628397E-4,
							-2.29549961613378126380E-4,
							-2.68132617805781232825E-3,
							3.47222221605458667310E-3,
							8.33333333333482257126E-2,
		};
            double MAXSTIR = 143.01608;

            double w = 1.0 / x;
            double y = Math.Exp(x);

            w = 1.0 + w * polevl(w, STIR, 4);

            if (x > MAXSTIR)
            {
                /* Avoid overflow in Math.Pow() */
                double v = Math.Pow(x, 0.5 * x - 0.25);
                y = v * (v / y);
            }
            else
            {
                y = Math.Pow(x, x - 0.5) / y;
            }
            y = SQTPI * y * w;
            return y;
        }

        /// <summary>
        /// Evaluates polynomial of degree N
        /// </summary>
        /// <param name="x"></param>
        /// <param name="coef"></param>
        /// <param name="N"></param>
        /// <returns></returns>
        private static double polevl(double x, double[] coef, int N)
        {
            double ans;

            ans = coef[0];

            for (int i = 1; i <= N; i++)
            {
                ans = ans * x + coef[i];
            }

            return ans;
        }

        private const double SQTPI = 2.50662827463100050242E0;
        private const double MACHEP = 1.11022302462515654042E-16;
        private const double MAXLOG = 7.09782712893383996732E2;
        private const double MINLOG = -7.451332191019412076235E2;
        private const double MAXGAM = 171.624376956302725;
        private const double LOGPI = 1.14472988584940017414;
    }
}
